---
title: "Figure 3"
output: 
---

# Figure 3

# Who's to Blame? Postconflict Violence and Public Attitudes Towards Peace Agreements
# Wyer, Frank. 

#clear environment
```{r clear environment}
rm(list = ls())
```

# uncomment and set working directory to replication archive
# setwd("~/blame_replication")

# Uncomment to install packages if necessary
# install.packages("tidyverse")
# install.packages(estimatr)

#load packages
```{r}
library(tidyverse)
library(estimatr)
```

#read in data
```{r}
survey_clean <- read.csv("survey_clean.csv")
```

#scale mechanism variables with respect to control group mean and sd
```{r scale mechanism variables}
#confidence in government implementation
survey_clean$govconf_scale <- (survey_clean$Q45 - mean(survey_clean$Q45[survey_clean$treatment == "C"], na.rm = TRUE)) / sd(survey_clean$Q45[survey_clean$treatment == "C"], na.rm = TRUE)

#confidence in rebel compliance
survey_clean$farcconf_scale <- (survey_clean$Q46 - mean(survey_clean$Q46[survey_clean$treatment == "C"], na.rm = TRUE)) / sd(survey_clean$Q46[survey_clean$treatment == "C"], na.rm = TRUE)
```

#estimate effects of treatments on mechanism questions
```{r estimate mechanism models}
implement_mechanism_95 <- lm_lin(formula = govconf_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)

implement_mechanism_90 <- lm_lin(formula = govconf_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff +  factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1) 

comply_mechanism_95 <- lm_lin(formula = farcconf_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .05)

comply_mechanism_90 <- lm_lin(formula = farcconf_scale ~ treatment, covariates = ~ Q15 + Q25 + urbandummy + engage_zscale + farc_presence + homratediff + factor(regionname), se_type = "HC2", data = survey_clean, alpha = .1)
```

#extract model coefficients to dataframe
```{r create dataframe of mechanism results with Confidence Gap}
mechanisms_df <- rbind(

data.frame(Treatment = "Government Culpability Treatment", Method = "Select Covariates", Conf_Level = "95%", Outcome = "Confidence in Government Implementation", Estimate = implement_mechanism_95$coefficients['treatmentT2A'], Conf_Low = implement_mechanism_95$conf.low['treatmentT2A'], Conf_High = implement_mechanism_95$conf.high['treatmentT2A']), 

data.frame(Treatment = "Government Culpability Treatment", Method = "Select Covariates", Conf_Level = "90%", Outcome = "Confidence in Government Implementation", Estimate = implement_mechanism_90$coefficients['treatmentT2A'], Conf_Low = implement_mechanism_90$conf.low['treatmentT2A'], Conf_High = implement_mechanism_90$conf.high['treatmentT2A']),   
    
data.frame(Treatment = "Rebel Culpability Treatment", Method = "Select Covariates", Conf_Level = "95%", Outcome = "Confidence in Government Implementation", Estimate = implement_mechanism_95$coefficients['treatmentT2B'], Conf_Low = implement_mechanism_95$conf.low['treatmentT2B'], Conf_High = implement_mechanism_95$conf.high['treatmentT2B']), 

data.frame(Treatment = "Rebel Culpability Treatment", Method = "Select Covariates", Conf_Level = "90%", Outcome = "Confidence in Government Implementation", Estimate = implement_mechanism_90$coefficients['treatmentT2B'], Conf_Low = implement_mechanism_90$conf.low['treatmentT2B'], Conf_High = implement_mechanism_90$conf.high['treatmentT2B']), 

data.frame(Treatment = "Government Culpability Treatment", Method = "Select Covariates", Conf_Level = "95%", Outcome = "Confidence in Rebel Compliance", Estimate = comply_mechanism_95$coefficients['treatmentT2A'], Conf_Low = comply_mechanism_95$conf.low['treatmentT2A'], Conf_High = comply_mechanism_95$conf.high['treatmentT2A']), 

data.frame(Treatment = "Government Culpability Treatment", Method = "Select Covariates", Conf_Level = "90%", Outcome = "Confidence in Rebel Compliance", Estimate = comply_mechanism_90$coefficients['treatmentT2A'], Conf_Low = comply_mechanism_90$conf.low['treatmentT2A'], Conf_High = comply_mechanism_90$conf.high['treatmentT2A']),   
    
data.frame(Treatment = "Rebel Culpability Treatment", Method = "Select Covariates", Conf_Level = "95%", Outcome = "Confidence in Rebel Compliance", Estimate = comply_mechanism_95$coefficients['treatmentT2B'], Conf_Low = comply_mechanism_95$conf.low['treatmentT2B'], Conf_High = comply_mechanism_95$conf.high['treatmentT2B']), 

data.frame(Treatment = "Rebel Culpability Treatment", Method = "Select Covariates", Conf_Level = "90%", Outcome = "Confidence in Rebel Compliance", Estimate = comply_mechanism_90$coefficients['treatmentT2B'], Conf_Low = comply_mechanism_90$conf.low['treatmentT2B'], Conf_High = comply_mechanism_90$conf.high['treatmentT2B'])
)
```

#convert to factors to control order of outcomes in plot
```{r control order of outcomes and treatments for mechanism plot}
mechanisms_df$Outcome <- factor(mechanisms_df$Outcome, levels = c("Confidence in Rebel Compliance", "Confidence in Government Implementation")) 

mechanisms_df$Treatment <- factor(mechanisms_df$Treatment, levels = c( "Government Culpability Treatment", "Rebel Culpability Treatment"))
```

#save mechanisms plot 
```{r mechanisms plot}
mech_plot_bw <- 
  ggplot(mechanisms_df %>% pivot_wider(names_from = "Conf_Level", values_from = c(Estimate, Conf_Low, Conf_High)), aes(y = Outcome, x = `Estimate_95%`, shape = Treatment)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +
  scale_shape_manual(name = "", values = c(15, 19), breaks = c("Rebel Culpability Treatment", "Government Culpability Treatment"), labels = function(x) str_wrap(x, width = 5)) +
  geom_linerange(aes(y = Outcome, xmin = `Conf_Low_90%`, xmax = `Conf_High_90%`), lwd = 1, position=position_dodge(.5)) + 
  geom_linerange(aes(y = Outcome, xmin = `Conf_Low_95%`, xmax = `Conf_High_95%`), lwd = 1/2, position=position_dodge(.5)) +
  geom_point(aes(y = Outcome, x = `Estimate_95%`), position=position_dodge(.5), size = 2.5) + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 15)) +
  guides(shape = guide_legend(byrow = TRUE)) + 
  labs(x = "", y = "", shape = "") +
  theme_bw() +
  theme(legend.key.size = unit(2.5, 'lines'))
```
